load("het_data.rdata")

## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(rdrobust)
library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)

mean_eff <- array(NA, dim = c(4, 2) )
mean_se  <- array(NA, dim = c(4, 2) )
mean_eff_prob <- array(NA, dim = c(4, 2) )
upp_prob  <- array(NA, dim = c(4, 2) )
low_prob  <- array(NA, dim = c(4, 2) )
years <- c(2009, 2013, 2014, 2015)


rownames(mean_eff) <- c("Child is in high school",
                        "Child is not in high school",
                        "Child is in high school rob",
                        "Child is not in high school rob")
colnames(mean_eff) <- c("Parent does have high school", "Parent does not have high school")

rownames(mean_se)  <- rownames(mean_eff)
colnames(mean_se)  <- colnames(mean_eff)

rownames(mean_eff_prob)  <- rownames(mean_eff)
colnames(mean_eff_prob)  <- colnames(mean_eff)
rownames(upp_prob)  <- rownames(mean_eff)
colnames(upp_prob)  <- colnames(mean_eff)
rownames(low_prob)  <- rownames(mean_eff)
colnames(low_prob)  <- colnames(mean_eff)

for (k in 1:2){
  for(j in 1:2){
    
    data <- 
      data.frame(data_list[[8]]) %>%
      filter(par_hs == k - 1 & child_hs == j - 1) %>%
      mutate(treated = three.days > 0)
    
    data_close <- 
      data %>%
      ungroup() %>%
      mutate(voted     = round(par_vote*n),
             abstained = round((1-par_vote)*n)) %>% 
      select(year, three.days, treated, voted, abstained)
    
    data_voted <- 
      as.data.frame(lapply(data_close, 
                           function(x,p) rep(x,p), 
                           data_close[["voted"]])) %>%
      mutate(voted = 1)
    
    data_abstained <- 
      as.data.frame(lapply(data_close, 
                           function(x,p) rep(x,p), 
                           data_close[["abstained"]])) %>%
      mutate(voted = 0)
    
    data_master <- 
      rbind(data_voted, data_abstained) %>% 
      filter(year != 2015)
    
    coef     <- rep(NA, 3)
    coef_rob <- rep(NA, 3)
    se       <- rep(NA, 3)
    se_rob   <- rep(NA, 3)
    
    #remove 2015, as educational data is not updated
    for (i in 1:3){
      sub_data <-
        data_master %>%
        filter(year == years[i])
      
      mod <- 
        with(sub_data, rdrobust(y = voted, x = three.days))
      
      coef[i]     <- mod$coef[1]
      coef_rob[i] <- mod$coef[3]
      se[i]       <- mod$se[1]
      se_rob[i]   <- mod$se[3]
    }
    
    
    mean_vote <- with(subset(data_master, three.days < 0), mean(voted))
    
    mean_eff[j,k] <- meta.summaries(coef, se)$summary
    mean_se[j,k]  <- meta.summaries(coef, se)$se.summary
    mean_eff_prob[j,k] <- 
      qnorm(mean_eff[j, k] + mean_vote) - 
      qnorm(mean_vote)
    upp_prob[j,k] <- 
      qnorm(mean_eff[j, k] + mean_vote + qnorm(0.975)*mean_se[j, k]) - 
      qnorm(mean_vote)
    low_prob[j,k] <- 
      qnorm(mean_eff[j, k] + mean_vote + qnorm(0.025)*mean_se[j, k]) - 
      qnorm(mean_vote)  
    
    mean_eff[j + 2,k] <- meta.summaries(coef_rob, se_rob)$summary
    mean_se[j + 2,k]  <- meta.summaries(coef_rob, se_rob)$se.summary
    mean_eff_prob[j + 2,k] <- 
      qnorm(mean_eff[j + 2,k] + mean_vote) - 
      qnorm(mean_vote)
    upp_prob[j + 2,k] <- 
      qnorm(mean_eff[j + 2,k] + mean_vote + qnorm(0.975)*mean_se[j + 2,k]) - 
      qnorm(mean_vote)
    low_prob[j + 2,k] <- 
      qnorm(mean_eff[j + 2,k] + mean_vote + qnorm(0.025)*mean_se[j + 2,k]) - 
      qnorm(mean_vote)  
    
    
    print(c(j, k)) 
  }
}


mean_eff <- mean_eff * 100
mean_se  <- mean_se  * 100

effs <- mean_eff[1,]
low  <- mean_eff[1,] + mean_se[1,] * qnorm(0.025) 
upp  <- mean_eff[1,] + mean_se[1,] * qnorm(0.975) 


for (i in 2:4){
  effs <- c(effs, mean_eff[i,])
  low  <- c(low, mean_eff[i,] + mean_se[i,] * qnorm(0.025) ) 
  upp  <- c(upp, mean_eff[i,] + mean_se[i,] * qnorm(0.975) )
  
}


for (i in 1:4){
  effs <- c(effs, mean_eff_prob[i,])
  low  <- c(low , low_prob[i,])
  upp  <- c(upp , upp_prob[i,])
}


labels <- c("Child is in high school  ",
            "Child is not in high school")

place  <- rep(c(4.5, 3.5, 2, 1), 4)
parent <- rep(c("Parent does have \nhigh school", "Parent does not \nhave high school"), 8)
type   <- rep(rep(c("Conventional RD", "Robust RD"), each = 4), 2)
prob   <- rep(c("Effect size in percentage points", "Effect size in probits"), each = 8)


plot_data <- 
  data_frame(effs, upp, low, place, parent, type, prob)


plot_prob <- 
  ggplot(data = plot_data,
         aes(y = place,
             x = effs,
             xmin = low,
             xmax = upp,
             fill = parent)) +
  geom_errorbarh(aes(alpha = parent, linetype = parent), height = 0) +
  facet_grid(type ~ prob, scales = "free_x") +
  scale_y_continuous(breaks = c(4, 1.5), labels = labels) +
  theme_classic() + 
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 16),
        legend.title= element_blank(), 
        legend.position = "top",
        panel.spacing = unit(2, "lines")) +
  ylab("") + 
  xlab("") +
  geom_point(aes(alpha = parent), size = 3) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  geom_vline(xintercept = 0, linetype = 3, alpha = 0.5)

# put baseline rates in table

tab <- matrix(NA, nrow = 4, ncol = 2)

#remove 2015, as educational data is not updated

for (i in 1:2){
  for (k in 1:2){
    data <- 
      data.frame(data_list[[8]]) %>%      
      filter(par_hs == i - 1 & year != 2015) %>%
      filter(three.days < 0)
    
    if(k == 1) data <- subset(data, child_hs == 0)
    if(k == 2) data <- subset(data, child_hs == 1)
    
    tab[7 - i - 2*k, 1] <- paste(round(weighted.mean(data$par_vote, data$n)*100, 2)) 
    tab[7 - i - 2*k, 2] <- paste("(",sum(data$n),")", sep = "")
  }
}

colnames(tab) <- c("Turnout", "N")
rownames(tab) <- c("Child in high school, parent with high school",
                   "Child in high school, parent without high school",
                   "Child not in high school, parent with high school",
                   "Child not in high school, parent without high school")

xtable(tab)


